home *** CD-ROM | disk | FTP | other *** search
/ NASA Climatology Interdisciplinary Data Collection / NASA Climatology Interdisciplinary Data Collection - Disc 4.iso / software / grads / doc / manual / corr.gs < prev    next >
Encoding:
Text File  |  1998-04-23  |  4.8 KB  |  256 lines
  1. function corr(args)
  2. dtg=subwrd(args,1)
  3. *
  4. * demonstration of calculating the correlation coef between
  5. * two horizontal (lon,lat) grids
  6. *
  7. * open a typical nwp file
  8. *
  9. if(dtg = '') ; 'quit' ; endif
  10. 'reinit'
  11. fname='/d1/obs/fnmoc/dat/nogaps.25.'dtg'.ctl'
  12. fn=ofile(fname)
  13. if(fn = 0 ) ; return ; endif
  14. *
  15. * set lat lon bounds
  16. *
  17. * NOTE: the underscore denotes "global" script variables
  18. * which do not have
  19. * to be passed to calling functions
  20. *
  21. _lat1=0
  22. _lat2=60
  23. _lon1=100
  24. _lon2=180
  25. *
  26. * set the dimension environment to 2-D lon,lat
  27. *
  28. 'set lat '_lat1' '_lat2
  29. 'set lon '_lon1' '_lon2
  30. 'set lev 500'
  31. 'set t 1'
  32. *
  33. * define script variables consisting of a string containing
  34. * the GrADS "expression" (a grid)
  35. *
  36. * in this case, compare the u comp at 850 with u at 200 mb
  37. *
  38. x1='u(lev=850)'
  39. x2='u(lev=200)'
  40. *
  41. * now try u at 500 between t=1 (analysis) and t=2 12 h forecast
  42. *
  43. *x1='u(t=1)'
  44. *x2='u(t=2)'
  45. *
  46. * calculate the unweighted and weighted corr coefficients
  47. *
  48. rhou=corru(x1,x2)
  49. rhow=corrw(x1,x2)
  50. *
  51. * contour plot of the fields; use q pos to continue after clicking
  52. * on the screen
  53. *
  54. 'c'
  55. 'set gxout contour'
  56. 'set ccolor rainbow'
  57. 'set cint 5'
  58. 'd 'x1
  59. 'draw title field #1'
  60. 'q pos'
  61. *
  62. * overlay second field
  63. *
  64. 'set ccolor 1'
  65. 'set cint 5'
  66. 'd 'x2
  67. 'draw title field #2'
  68. 'q pos'
  69. *
  70. * clear and look at a scatter plot
  71. *
  72. 'c'
  73. 'set gxout scatter'
  74. 'd 'x1';'x2
  75. *
  76. * put the coef in the scatter plot title, note the fancy use
  77. * of fonts and subscripts...
  78. *
  79. 'draw title `3r`0`barea weighted`n = 'rhow' `3r`0`bunweighted`n = 'rhou
  80. return
  81. function corrw(x1,x2)
  82. *
  83. * calc correlation coefficient for 2-D horizontal fields using aave
  84. *
  85. * for non horizontal fields use ave(ave...)
  86. *
  87. * ASSUMES THE GRIDS HAVE THE SAME NUMBER OF DEFINED POINTS!!!!
  88. *
  89. * based on formula on p.92 of Panofsky and Brier, 1968 "Some applications of
  90. * statistics to meteorology." Penn State Press
  91. *
  92. *
  93. * set dimension environment to 0-D (a point), as a consequence
  94. * all calcs are written to the string variable result
  95. *
  96. *'set x 1'
  97. *'set y 1'
  98. *
  99. * display (calc) the ave of field #1
  100. *
  101. 'd aave('x1',lon='_lon1',lon='_lon2',lat='_lat1',lat='_lat2')'
  102. ave1=subwrd(result,4)
  103. *
  104. * field #2
  105. *
  106. 'd aave('x2',lon='_lon1',lon='_lon2',lat='_lat1',lat='_lat2')'
  107. ave2=subwrd(result,4)
  108. *
  109. * cross product #1 * #2
  110. *
  111. 'd aave('x1'*'x2',lon='_lon1',lon='_lon2',lat='_lat1',lat='_lat2')'
  112. ave12=subwrd(result,4)
  113. *
  114. * ave of the sqr #1
  115. *
  116. 'd aave('x1'*'x1',lon='_lon1',lon='_lon2',lat='_lat1',lat='_lat2')'
  117. ave1s=subwrd(result,4)
  118. *
  119. * ave of sqr of #2
  120. 'd aave('x2'*'x2',lon='_lon1',lon='_lon2',lat='_lat1',lat='_lat2')'
  121. ave2s=subwrd(result,4)
  122. *
  123. * the first part of the denominator
  124. *
  125. 'd sqrt('ave1s'-'ave1'*'ave1')'
  126. den1=subwrd(result,4)
  127. *
  128. * the second part
  129. *
  130. 'd sqrt('ave2s'-'ave2'*'ave2')'
  131. den2=subwrd(result,4)
  132. *
  133. * the numerator
  134. *
  135. 'd 'ave12'-'ave1'*'ave2
  136. num=subwrd(result,4)
  137. *
  138. * the coefficient
  139. *
  140. 'd 'num'/('den1'*'den2')'
  141. corr=subwrd(result,4)
  142. *
  143. * output to the terminal
  144. *
  145. say ' The area weighted correlation = 'corr
  146. *
  147. * return the string to the main program
  148. *
  149. return(corr)
  150. function corru(x1,x2)
  151. *
  152. * use the new version of gxout stat to get the # of valid points
  153. * (for info purposes only)
  154. * and the biased (/N) vice the unbiased (/(N-1)) averages
  155. *
  156. * ASSUMES THE GRIDS HAVE THE SAME NUMBER OF DEFINED POINTS!!!!
  157. *
  158. * imax is the number of lines produced by gxout stat so the loop
  159. * doesn't go forever...
  160. *
  161. * based on formula on p.92 of Panofsky and Brier, 1968 "Some applications of
  162. * statistics to meteorology." Penn State Press
  163. *
  164. imax=12
  165. 'set gxout stat'
  166. 'd 'x1
  167. result1=result
  168. 'd 'x2
  169. result2=result
  170. 'd 'x1'*'x2
  171. result3=result
  172. *
  173. * parse the strings for the info
  174. *
  175. i=1
  176. while(i<=imax)
  177. card1=sublin(result1,i)
  178. card2=sublin(result2,i)
  179. card3=sublin(result3,i)
  180. *
  181. * get number of defined points
  182. *
  183. if(subwrd(card1,2)='count')
  184. npt=subwrd(card1,8)
  185. endif
  186. *
  187. * the sum and sumsqr / n (biased)
  188. *
  189. if(subwrd(card1,1)='Stats(sum,sumsqr)/n:')
  190. ave1b=subwrd(card1,2)
  191. avesq1b=subwrd(card1,3)
  192. ave2b=subwrd(card2,2)
  193. avesq2b=subwrd(card2,3)
  194. ave12b=subwrd(card3,2)
  195. avesq12b=subwrd(card3,3)
  196. *
  197. * break the loop now that we have what we need
  198. *
  199. break
  200. endif
  201. i=i+1
  202. endwhile
  203. *
  204. * turn off stat output and do the calculation using display
  205. *
  206. * the result goes to the script variable result
  207. *
  208. 'set gxout contour'
  209. 'd 'ave12b'-'ave1b'*'ave2b
  210. num=subwrd(result,4)
  211. 'd sqrt('avesq1b'-'ave1b'*'ave1b')'
  212. den1=subwrd(result,4)
  213. 'd sqrt('avesq2b'-'ave2b'*'ave2b')'
  214. den2=subwrd(result,4)
  215. 'd 'num'/('den1'*'den2')'
  216. corr=subwrd(result,4)
  217. say ' The UNWEIGHTED correlation = 'corr
  218. *
  219. * return to main program
  220. *
  221. return(corr)
  222. *
  223. *-------------------------- ofile ------------------
  224. *
  225. function ofile (fname)
  226. 'query files'
  227. i = 0
  228. while (1)
  229. if (subwrd(result,1)='No')
  230. ret = 0
  231. break;
  232. endif
  233. rec = sublin(result,i*3+2)
  234. if (rec='')
  235. ret = 0;
  236. break;
  237. endif
  238. if (subwrd(rec,2)=fname)
  239. rec = sublin(result,i*3+1)
  240. ret = subwrd(rec,2)
  241. break;
  242. endif
  243. i = i + 1
  244. endwhile
  245. if (ret=0)
  246. 'open 'fname
  247. if (rc>0)
  248. say "Error opening "fname
  249. return (0)
  250. endif
  251. rec = sublin(result,2)
  252. ret = subwrd(rec,8)
  253. endif
  254. return (ret)
  255.  
  256.